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SURFACE  RECONSTRUCTION  PRESERVING  DISCONTINUITIES. 


J.  L.  Marroquin 


ABSTRACT:  This  paper  presents  some  experimental  results  that  indicate 
the  plausibility  of  using  non-convex  variational  principles  to  reconstruct 
piecewise  smooth  surfaces  from  sparse  and  noisy  data.  This  method 
uses  prior  generic  knowledge  about  the  geometry  of  the  discontinuities  to 
prevent  the  blurring  of  the  boundaries  between  continuous  subregions. 

We  include  examples  of  the  application  of  this  approach  to  the  reconstruc¬ 
tion  of  synthetic  surfaces,  and  to  the  interpolation  of  disparity  data  from 
the  stereo  processing  of  real  images. 
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I.  Introduction. 


(n  ihc  processing  of  two-dimensional  signals,  often  one  lias  the  problem  of 
reconstructing  a  piecewise  smooth  surface  from  noisy  observations  taken  at  sparse 
locations.  '  In  this  reconstruction,  it  is  important  not  only  to  interpolate  smooth 
patches  over  uniform  regions,  but  to  locate  and  preserve  the  discontinuities  that 
bound  these  regions,  since  very  often  they  are  the  most  important  parts  of  the 
surface.  They  may  represent  object  boundaries  in  vision  problems  (such  as  image 
segmentation;  depth  from  stereo;  shape  from  shading;  structure  from  motion, 
etc.);  geological  faults  in  geophysical  information  processing,  etc. 

The  most  successful  approaches  to  this  problem  [T2]  consist  of,  first,  inter¬ 
polating  an  everywhere  smooth  surface  over  the  whole  domain;  then,  applying 
some  kind  of  discontinuity  detector  (followed  by  a  thresholding  operation)  to  try 
to  find  the  significant  boundaries,  and  finally,  to  re- interpolate  smooth  patches 
over  the  continuous  subregions. 

The  results  that  have  been  obtained  with  this  technique,  however,  are  not 
completely  satisfactory.  The  main  probtem  is  that  the  task  of  the  discontinuity 
detector  is  hindered  by  the  previous  smooth  interpolation  operation.  This  becomes 
critical  when  the  observations  are  sparsely  located,  since  in  this  case,  the  discon¬ 
tinuities  may  be  smeared  in  the  interpolation  phase  to  such  a  degree  that  it  may 
become  impossible  to  recover  them  in  the  detection  phase. 

One  way  around  this  difficulty  is  to  perform  the  boundary  detection  and 
interpolation  tasks  at  the  same  time.  In  the  method  we  will  present,  this  is  done 
by  generating  a  variational  principle  that  includes  our  prior  knowledge  about  the 
smoothness  of  the  surface  and  about  the  geometry  of  the  discontinuities,  as  well 
as  the  information  provided  by  the  observations.  The  global  minimum  of  this 
(non-convex)  "energy"  functional  is  then  found  by  a  stochastic  approximation 
scheme. 

This  approach  is  based  on  the  work  of  the  Geman  brothers  [01].  Before 
describing  it,  let  us  formulate  the  problem  in  a  more  precise  way. 

Imagine  a  region  n  of  the  plane  which  is  formed  by  a  number  of  subregions 
separated  by  boundaries  which  are  known  to  be  piecewise  smooth  curves.  Suppose 
that  within  each  of  these  subregions,  some  property  /  (in  what  follows,  we  will  refer 
to  /  as  "depth")  varies  in  a  smooth  fashion,  presenting,  at  the  same  time,  abrupt 
jumps  across  most  of  the  boundaries.  Suppose  also  that  we  have  measurements 
for  the  values  of  /  at  some  discrete  set  of  sites  5;  these  measurements  will,  in 
general,  be  corrupted  by  some  form  of  noise. 

Our  problem  is  then  to  estimate  the  values  of  /  on  some  finite  lattice  of 
points  L  C  Cl,  and  to  find  the  position  of  the  boundaries,  using  all  the  available 
information  in  an  optimal  way. 
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2.  Genian's  Work  on  Bayesian  Image  Restoration. 

2.1  Images  as  Markov  Random  Fields. 

The  first  idea  on  which  this  approach  is  based,  is  that  an  image  formed 
by  regions  of  constant  intensity  separated  by  piecewise  smooth  boundaries  can 
be  modeled  as  a  sample  function  of  a  stochastic  process  based  on  the  Gibbs 
distribution,  that  is,  as  a  Markov  Random  Field  (MRF)  (sample  functions  of 
actual  MRF's  can  be  found  in  his  paper.  See  also  [Cl]  and  [G2]).  This  means 
that  the  conditional  probability  of  a  given  pixel  j  having  a  particular  value  /y, 
given  the  values  of  /  in  all  the  remaining  sites  of  the  lattice,  is  identical  to  the 
conditional  probability  given  the  values  of  /  in  a  small  set  of  sites  which  we  will 
call  the  neighborhood  of  j. 

Given  a  system  of  neighborhoods  on  a  lattice,  we  define  a  "clique"  C  as  a 
set  of  sites  such  that  all  the  sites  that  belong  to  C  are  neighbours  of  each  other. 
For  example,  on  a  4-connccted  lattice  (Fig.  1-a),  the  sites  1,  2,  3  and  4  form  the 
neighborhood  of  site  j ,  and  the  cliques  are  sets  consisting  either  of  single  sites,  or 
of  two  (vertically  or  horizontally)  adjacent  sites  (nearest  neighbours;  see  fig.  1-b), 

It  can  be  shown  [PI],  that  if  /  is  a  MRF  with  a  given  system  of  neighborhoods, 
the  probability  density  of  the  images  generated  bv  this  process  is  nf  the  form: 

pad = ?(“  =  /) = 

where  Z  is  a  normalizing  constant,  /?  is  a  parameter,  and  the  "Energy  function" 
U{f)  is  of  the  form; 

tf(/)=x;w) 

c 

where  the  potentials  Vc(f)  are  functions  supported  on  the  cliques  associated  to 
the  given  neighborhood  system.  Thus,  in  our  example  of  a  4-connected  lattice,  U 
would  be: 

rt/)-  EWi)+  £  vHfi.fi) 

j  i,]ENN 

where  i,j  6  Ns  means  that  i  and  j  are  nearest  neighbours,  and  Vx  and  V2  are 
some  functions. 

In  particular,  if  we  want  a  Markov  random  process  that  generates  piecewise 
constant  surfaces,  we  may  use  potentials: 

Vi(/)  =  0 
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Fig.  1-a.  Sites  1,  2,  3  and  4  are  the  neighborhood  of  site  j. 
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Fig.  1-b.  Cliques  for  the  4-connected  lattice  of  Fig.  1-a 
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The  parameter  0  can  be  interpreted  as  the  natural  temperature  of  the  system, 
and  controls  the  expected  size  of  the  regions  formed  by  the  process,  larger  regions 
being  formed  at  low  temperatures. 

2.2  Bayesian  Estimation  of  Markov  Random  Melds. 


Now,  suppose  that  we  have  observations  g  that  can  be  modeled  as: 

9j  =  ¥(//,■(/),  «y)  j€S  (2) 

where  ny  is  a  white  noise  process  uncorrelated  with  /;  S  is  some  set  of  sites;  Hj 
is  an  operator  with  local  support  (representing  blurring,  for  example),  and  is 
invertible  with  respect  to  ny.  (W  may  represent,  for  example,  noise  addition  or 
multiplication  followed  by  a  memoryless  transformation). 

The  conditional  probability  P(g  |  /)  is  given  by: 

P(9\f)=  II  Pn(*-\gj  -  Hj(f)) 

j€S 

where  Pn  is  the  probability  density  function  of  the  noise. 

The  posterior  distribution  is  found  by  Bayes  rule: 

p,u)  ■  ni  i  /) 


pu  i  ?)  = 


p(i) 


Replacing  the  expressions  for  P/(f)  and  P{o  |  /),  taking  logarithms,  and 
remembering  that  P(ff)  is  a  constant  for  a  given  set  of  observations,  we  get  that 
the  Maximum  a  Posteriori  (MAP)  estimate  for  /  is  found  by  minimizing: 


£(/)  -  \v(t)  -  £  l»[P.(*-1(w  -  //,(/))! 
p  yes 


(3) 


In  particular,  if  n  is  a  zero-mean  Gaussian  white  noise  stationary  process  with 
power  spectral  density  a2,  and 


9j  =  H,{f)  +  ny  =  fj  +  ny  , 


then. 
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P(3  I  /)  =  ___exp[___  22 (gj  -  fj)2] 


\f2io  60  '  j€S 

For  our  example  of  piecewise  constant  surfaces,  with  potentials  given  by  (1), 
the  MAP  estimate  will  be  obtained  by  minimizing: 


£(/)=  £  + A  £(/,-*,■)“ 

ij6Ns  la  jeS 


(4) 


As  we  can  see,  this  expression  has  two  terms:  One  that  measures  the 
agreement  of  the  estimate  with  the  observations,  and  another  that  corresponds  to 
the  constraints  imposed  by  our  prior  knowledge  about  the  nature  of  the  solution. 
The  tradeoff  between  them  is  controlled  ''y  tb"  parameter 

0 

“  2a2 

which  corresponds  to  the  signal  to  noise  ratio. 

It  is  interesting  to  note  that  the  general  form  of  this  expression  is  similar  to 
the  variational  principles  obtained  by  regularization  methods  for  solving  ill-posed 
problems,  a  playing  the  role  of  the  regularization  parameter.  In  the  standard 
regularization  methods,  however,  the  functional  is  convex  (because  of  the  choices 
of  norms  and  stabilizing  functionals,  see  [P2,T1]),  whereas  in  the  case  of  equation 
(4),  the  high  non-linearity  of  V2  makes  E(f)  non-convex.  For  this  reason,  its 
minimization  becomes  computationally  much  more  expensive. 

2.3  Simulated  Annealing  and  the  Minimization  of  E(f). 

In  a  recent  paper,  Kirkpatrick,  et.  al.  [Kl]  proposed  a  stochastic  approximation 
method  for  solving  combinatorial  optimization  problems,  which  can  be  used  for 
our  minimization  problem.  It  is  based  on  an  algorithm  invented  by  Metropolis 
[Ml]  that  simulates  the  behaviour  of  many-particle  systems  in  thermal  equilibrium: 

Consider  a  system  with  N  particles,  each  of  which  may  be  in  any  one  of 
a  discrete  number  of  allowable  states.  Let  /y  denote  the  current  state  of  the 
jth  particle;  T,  the  temperature,  and  let  E(f)  be  the  total  energy  of  the  system. 
Suppose  that  we  visit  the  particles  of  the  system  in  some  random  sequential  order. 
When  a  particle  j  is  visited,  we  update  its  state  as  follows: 

(i)  Choose  a  new  state  /y  randomly  from  the  set  of  allowable  states  (excluding 
the  current  state  /,) ,  using  a  uniformly  distributed  random  number. 

(ii)  Compute  the  increment  in  energy  A£y  that  results  from  moving  the  state 
of  the  jth  particle  from  fj  to  fj. 

(iii)  If  A25y  <  0,  make  the  move,  i.e.,  set  fj  —  fj. 

If  AEj  >  0,  generate  a  new  random  number  r,  uniformly  distributed 
between  0  and  1. 

If  r  <  e~±E*/T,  set  fj  =  /y. 

If  r  >  e_A®»/r,  leave  /y  unchanged. 

It  can  be  shown  [M1,K2]  that  if  every  particle  is  visited  infinitely  often, 
this  procedure  will  generate  global  states  for  /  distributed  according  to  Gibbs 
distribution,  i.e.,  as  the  number  of  iterations  goes  to  infinity,  we  will  have: 
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As  T  goes  to  0,  this  distribution  will  tend  to  an  impulse  (or  set  of  impulses) 
corresponding  to  the  state  (or  states)  of  minimum  energy,  that  is,  to  the  value  of 
/  that  minimizes  E(f)  globally. 

One  serious  difficulty,  however,  is  that  attaining  thermal  equilibrium  might 
take  a  very  long  time  at  low  temperatures.  Kirkpatrick's  idea  was  to  start  at 
a  relatively  high  temperature  (where  thermal  equilibrium  is  reached  very  fast), 
and  then,  to  slowly  cool  the  system,  until  "freezing”  occurs  and  the  state  stops 
changing. 

Geman  &  Geman  were  able  to  show  that  if  the  temperature  is  lowered  at  the 

rate: 

T  =  log(n  +  1)  (5) 

where  n  is  the  number  of  iterations,  and  C  is  a  constant,  this  algorithm  will  in 
fact  converge  (in  probability)  to  the  set  of  states  of  minimal  energy  (see  also  (G3J). 
Also,  they  proved  that  the  non-homogeneous  Markov  chain  that  corresponds  to 
the  annealing  procedure  is  (asymptotically)  stationary  and  ergodic,  so  that  we  may 
use  time  statistics  to  estimate  the  final  state.  We  will  return  to  this  in  the  next 
section. 

Unfortunately,  the  value  of  the  constant  C  necessary  to  guarantee  convergence 
is  in  general  very  high  (so  that  the  convergence  time  becomes  i m p ractica  1 1  y  slow), 
but  it  has  been  found  experimentally  that  a  value  of  C  —  3/3  log  2  (where  (3  is  the 
natural  temperature  of  the  system)  normally  produces  a  reasonable  convergence 
behaviour  (of  the  order  of  1000  iterations). 

The  computational  efficiency  of  Kirkpatrick’s  algorithm  depends  on  whether 
the  increment  in  energy  A  Ej  associated  to  a  change  in  the  state  of  the  jUl  variable 
is  easy  to  compute.  Fortunately,  the  Markov  property  of  /,  and  the  localization 
condition  on  the  support  of  the  operator  H}  (equation  (2))  guarantee  that,  for  our 
problem,  this  will  be  the  case,  so  that  simulated  annealing  is  a  practical,  although 
expensive  way  of  finding  the  minimum  of  (3). 


2.4  The  Line  Process. 


Another  powerful  idea  in  the  Geman’s  work  is  the  introduction  of  an 
unobservable  line  process  l  into  the  image  model.  The  variables  associated  with 
this  process  are  located  at  the  sites  of  the  dual  lattice  of  lines  that  connect  the 
sites  of  the  original  (pixel)  lattice  (see  Fig.  2-a).  These  variables  may  be  binary 
(indicating  the  presence  or  absence  of  a  boundary  between  two  pixels),  or  may 
take  more  values  to  indicate  the  orientation  of  the  boundary  as  well.  In  both 
cases,  their  function  is  to  decouple  adjacent  pixels,  reducing  the  total  energy  if  the 
intensities  of  these  pixels  are  different. 


OXOXOXOXOXOXO 

x  x  x  x  x  x  x 
OXOXOXOXOXOXO 
x  x  x  x  x  x  x 
OXOXOXOXOXOXO 
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OXOXOXOXOXOXO 
x  x  x  x  x  x  x 
OXOXOXOXOXOXO 


Fig.  2-a.  Dual  Lattice  of  Line  elements  (sites  denoted  by  x). 
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Fig  2-b.  Giques  for  the  line  process. 
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The  introduction  of  this  process  is  particularly  important  for  the 
discontinuity-preserving  reconstruction  of  surfaces  from  sparse  observa¬ 
tions,  since  it  allows  us  to  include  into  the  estimation  problem  our  prior 
knowledge  about  the  geometry  of  the  boundaries  between  regions  in  an 
explicit  way.  This  is  done  by  modifying  the  energy  function;  the  new 
expression  is: 


£(/)=  £  +  + 


C, 


jes 


(6) 


where 


V,(f  f-  l-  )  =  J®’  **  IS  on 

tj)  (/,,/>),  otherwise 


if  lij  is 


Vo  is  defined  in  eq.  (1);  ltj  is  the  line  element  between  sites  i  and  j, 
and  the  line  potentials  Vc,  have  as  supports  cliques  of  size  4,  such  as  the 
one  shown  in  Fig.  2-b.  Every  line  element  (except  at  the  boundaries 
of  the  lattice)  belongs  to  2  such  cliques.  The  values  of  the  potentials 
associated  with  each  possible  configuration  of  lines  within  a  clique  must 
be  specified.  Thus,  for  example,  if  we  know  that  straight  horizontal  and 
vertical  boundaries  are  likely  to  be  present,  we  may  use  a  binary  process, 
and  potential  values  as  those  of  Fig.  3  (rotational  invariance  is  assumed). 
If  we  want  to  handle  more  general  situations  (such  as  piecewise  smooth 
boundaries),  it  is  necessary  to  allow  more  states  for  the  line  elements, 
corresponding  to  different  orientations,  augmenting  consequently  the  table 
of  values  for  tire  potentials. 


3.  Extension  to  the  Boundary  Preserving  Interpolation  of  Piecewise 
Smooth  Surfaces. 


3.1  Continuous  Intensity  Values. 

Geman  &  Geman  successfully  applied  their  method  to  the  restoration 
of  piecewise  constant  images  corrupted  by  additive  white  Gaussian  noise, 
assuming  that  the  grey  levels  corresponding  to  the  constant  subregions 
were  either  known  in  advance,  or  easily  obtainable  from  the  data  (as 
peaks  of  a  histogram,  for  example).  However,  if  we  want  to  apply  this 
construction  to  more  practical  cases,  we  must  relax  the  piecewise  constant 
condition,  so  as  to  include  tilted  planes,  smooth  gradients,  and  in  general, 
piecewise  smooth  surfaces. 

To  do  this,  we  must  allow  the  depth  variables  fj  to  take  any  real 
value.  However,  this  would  require,  in  general,  to  replace  Metropolis 
algorithm  by  a  process  capable  of  handling  continuous  variables  (such 


as  a  diffusion.  See  [G2]),  and  in  this  case,  the  convergence  of  the 
annealing  process  is  uncertain  (rigorous  convergence  results  have  not 
been  established  yet). 

One  way  around  this  difficulty,  is  to  define  a  "mixed"  annealing 
strategy,  by  replacing  the  depth  potential  Vj  in  eq.  (6)  by  one  that 
guarantees  that,  for  any  given  state  of  the  line  process,  the  resulting 
conditional  energy  function  E(f  |  l)  is  convex,  so  that  the  use  of  a 
deterministic  gradient  descent  iteration  for  updating  the  state  of  the  depth 
process  is  justified  in  some  sense  (this  is  not  a  rigorous  argument;  we 
present  it  only  to  provide  some  motivation  for  the  definition  of  this 
strategy  which,  at  this  point,  is  justified  only  by  our  experimental  results). 

This  can  be  achieved  by  using  a  quadratic  potential: 


V/ifiJjJij 


if  kj  is  "on" 
otherwise 


(7) 


Since  the  resulting  conditional  energy  function  is  of  the  form: 

bu  1 1)  =  E( A  -  fi?  +  A  E  (/>  -  mf  +  K  («) 

la  jes 

for  some  positive  constant  K,  we  have  that  E{f  |  l)  >  0  for  any 
/  7^  0,  and  by  an  appropriate  change  of  coordinates,  it  can  be  put  in  the 
form: 

E{f  |  I)  =  utAu 

with  u  €  R},jK  and  A  a  non-negative  definite  matrix. 

Now,  for  any  t  e  (0, 1),  and  any  u^«,  we  have, 

tE{u  1 1)  +  (1  -  1 1)  -  E(tu  +  (1  -  t)v  |  i)  = 

=  t(l  —  t)(u  —  v)tA[u  —  v)  >  0 

so  that  E(f  j  ()  is  a  convex  function  of  /,  and  an  iteration  of  a  gradient 
descent  algorithm  will  move  towards  a  global  minimum  of  this  function. 

Note  however,  that  there  may  be  degenerate  cases  in  which  some  region 
Q  within  which  there  are  no  observations  is  isolated  from  the  rest  of  the 
lattice  by  the  line  process.  In  this  case,  any  solution  for  which 

fj  =  constant ,  ;  G  Q 


will  be  a  fixed  point  of  the  gradient  descent  algorithm,  and  although  all 
these  solutions  have  the  same  energy,  some  of  them  are  obviously  more 
desirable  than  others.  Experimentally,  we  have  found  that  a  good  strategy 


to  prevent  the  formation  of  undesirable  "islands”,  is  to  use  the  global 
minimum  of  E{f  |  /  =  0)  (which  is  the  unique  fixed  point  of  Lhe  gradient 
descent  algorithm)  as  the  starting  configuration  for  the  mixed  annealing 
process.  This  configuration  can  be  interpret'd  physically  as  a  membrane 
that  is  coupled  to  the  data  points  (observations)  by  means  of  springs  with 
constants  equal  to  r^,  and  is  in  a  position  that  minimizes  the  total  energy 
stored  in  itself,  and  in  the  springs  (see  f T2]). 

3.2  Tlie  Mixed  Annealing  Process. 

The  scheme  we  are  proposing  is  as  follows: 

Every  global  iteration  (that  is,  for  every  fixed  temperature),  all  the  line 
and  intensity  sites  are  visited  sequentially.  When  a  line  site  is  visited,  its 
state  is  updated  using  Metropolis  algorithm.  The  corresponding  increment 
in  energy  AE{  is  computed  as  follows: 

Let  Ci  and  be  the  two  cliques  to  which  the  line  element  belongs. 

Let  lij  be  its  current  state,  and  l,y  be  the  candidate  state  (if  /  is  a  binary 
process,  =  1  —  lij",  otherwise,  it  is  chosen  at  random  from  the  set 
of  allowable  states  different  from  /,y).  Let  l  be  the  current  global  line 
configuration,  and  /  be  obtained  from  l  by  replacing  ltJ  by  /tJ.  We  have: 

AE,  =  VCl(l)  ~  VCl( 0  +  VcX)  ~  VCl{ l)  + 

+Sgn(/,y-],y)(/t-/y)2  (9) 

where 

fi.  if  x  >  0 

sgn(i)  =  j  -l,  if  x  <  0 
lo,  if  x  =  0 

When  an  intensity  site  j  is  visited,  its  new  state  /’y  is  obtained 
deterministically  by  the  formula: 

_  ^»=l»-j|€(0,ll(1  ~  s6n(^y))/»  +  a<Jj9j 
£«:|»'-y|e(0,i|(l  —  Sgn(/,y))  +  Qr«7y 

where  a  =  and  qy  =  0,  unless  there  is  an  observation  at  site  j, 
in  which  case,  qy  =  1. 

The  temperature  is  lowered  using  (5). 

3.3  Parameter  Values. 

Unfortunately,  we  do  not  have,  at  this  point,  a  clean  way  for  selecting 
the  values  for  the  parameters  of  the  system,  other  than  a  trial  and  error 
procedure.  A  few  considerations  about  their  meaning  are  in  order: 
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The  parameter  a  controls  the  degree  of  smoothing  that  one  applies 
to  the  data,  and  must  be  related  to  the  quality  of  the  observations.  Large 
values  of  a  («=»  10)  will  force  the  interpolated  surface  to  pass  through 
the  observations,  while  small  values  (^  l)  will  smooth  away  all  but  the 
largest  fluctuations. 

The  parameters  associated  with  the  line  potentials  control  both  the 
shape  and  the  number  of  boundaries  that  the  algorithm  will  detect.  For 
example,  for  the  binary  line  process  of  Fig.  3,  if  h  is  the  height  of  the 
smallest  jump  that  we  want  to  consider  as  a  boundary,  and  d  is  the  largest 
gradient  in  the  smooth  regions,  we  must  have  for  V\  (the  potential  of  a 
straight  boundary) 

d?  <V i  <  h 2  (11) 

The  ratio  of  V\  to  the  potentials  for  the  rest  of  the  configurations, 
represents  our  prior  knowledge  about  the  relative  likelihood  of  corners, 

"T”  junctions,  etc. 

4.  Experiments 

We  now  present  some  experimental  results  that  support  the  use  of 
this  approach  for  surface  reconstruction  and  image  segmentation  Lusks.  In 
these  examples,  we  assume  that  we  have  the  following  prior  know  ledge 
about  the  nature  of  the  surfaces  we  are  trying  to  reconstruct: 

(i)  The  region  under  consideration  can  be  segmented  into  a  small 
number  of  subregions. 

(ii)  Within  each  subregion  the  surface  is  smooth  (the  gradient  is  less 
than  0.5). 

(iii)  The  boundaries  between  regions  are  piecewise  horizontal  or 
vertical.  There  are  relatively  few  corners. 

(iv)  The  average  height  of  the  discontinuities  across  boundaries  is 
greater  than  0.8. 

(v)  The  observations  are  corrupted  by  an  additive  white  Gaussian 
noise  process,  and  we  have  some  estimate  of  its  intensity. 

This  knowledge  is  embodied  in  the  model  for  the  line  process,  and 
in  the  numerical  value  of  the  parameters.  For  our  experiments,  we  have 
used  the  binary  model  of  Fig.  3,  with  the  values  for  the  potentials  for 
the  different  configurations  indicated  there. 

4.1  Experimental  Results. 

In  the  first  set  of  experiments,  wc  generated  sparse  observation  points 
at  200  random  locations  of  a  30  X  30  rectangular  grid.  Figures  4,  5,  6 
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and  7  show  (with  height  coded  by  grey  level)  die  observations  (a);  the 
boundaries  found  by  the  algorithm  (b);  the  configuration  obtained  by 
interpolation  with  no  boundaries  (c),  and  the  final  reconstructed  surface 
(d),  for: 

(i)  A  square  at  height  2.0  over  a  background  at  constant  height  = 

1.0  (Fig.  4). 

(ii)  A  triangle,  with  the  same  characteristics  (Fig.  5). 

(iii)  A  tilted  square  plane  (slope  =  0.1)  over  a  constant  height 
background  with  white  Gaussian  added  noise  {a  =  O.lXFig.  6). 

(iv)  Three  rectangles  at  different  (constant)  heights  over  a  uniform 
background  (Fig.  7). 


Fig.  4.  (a)  Observations  of  a  square  at  height  2.0  over  a  back¬ 
ground  at  height  1.0.  (b)  Boundaries  found  by  the  Algorithm, 
(c)  Interpolation  with  no  boundaries,  (d)  Reconstructed  surface. 


Fig.  6.  (a)  Observations  of  a  tilted  square  (slope  =  0.1)  over  a 
background  at  height  1.0  with  added  white  Gaussian  noise  (a  = 
0.1)  (b)  Boundaries  found  by  the  Algorithm,  (c)  Interpolation 
with  no  boundaries,  (d)  Reconstructed  surface. 


Fig.  7.  (a)  Observations  of  3  rectangles  at  heights  2.0,  2.0 
and  3.0  over  a  background  at  height  1.0.  (b)  Boundaries  found 
by  the  Algorithm,  (c)  Interpolation  with  no  boundaries,  (d) 
Reconstructed  surface. 
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In  the  second  set  of  experiments,  the  same  algorithm  was  used  for 
a  boundary  detection  /  image  segmentation  task.  In  this  case,  we  have 
observations  (corrupted  by  white  Gaussian  noise)  at  almost  every  point 
in  the  lattice.  The  original  figure  is  a  10  X  10  tilted  square  plane  of  slope 
0.2  located  at  the  center  of  the  lattice  (Fig.  8).  Note  that  in  this  case,  the 
tilted  plane  cuts  across  the  uniform  background,  so  that  the  vertical  steps 
at  both  sides  have  opposite  signs,  while  the  horizontal  steps  change  sign 
at  the  center  of  the  figure,  where,  in  fact,  there  is  no  discontinuity.  As  in 
the  previous  experiments,  the  results  speak  for  themselves. 

Finally,  we  present  an  example  of  the  application  of  this  algorithm 
to  the  processing  of  real  images.  We  use  it  to  interpolate  disparity  data, 
obtained  along  the  zero-crossing  contours  of  the  convolution  of  a  stereo 
pair  of  aerial  photographs  with  a  "Difference  of  Gaussians"  operator,  by 
Grimson’s  implementation  of  the  Marr-Poggio  stereo  algorithm  [G4,M2] 
(the  data  array  was  rotated,  so  that  assumption  (iii)  held).  The  results  are 
shown  in  figure  9.  We  believe  that  they  will  improve  when  we  implement 
the  extensions  to  this  method  outlined  in  section  5.2. 

4.2  Convergence  of  the  Annealing  Process. 

Although  the  mixed  annealing  algorithm,  with  annealing  schedule 
given  by  (5),  eventually  converges  to  a  low  energy  solution,  we  have 
found  that  convergence  can  be  greatly  improved  if  we  periodically  set 
the  state  of  the  line  process  using  an  estimate  of  the  final  (lowest  energy) 
state. 

To  get  this  estimate,  we  use  the  ergodicity  of  the  process  ([G1J, 
theorem  C),  and  compute  time  statistics  about  its  evolution.  Specifically, 
we  estimate  the  marginal  distributions  of  the  states  of  every  "line"  cell  by 
computing,  within  a  time  window,  the  percentage  of  the  time  for  which 
this  cell  is  "on".  If  it  is  more  than  half  of  the  time,  the  state  of  this  cell 
in  the  annealing  network  is  forced  to  1;  otherwise,  it  is  forced  to  0.  As 
the  temperature  decreases,  the  probability  distribution  of  the  global  states 
becomes  peaked,  and  the  maximum  probability  estimate  obtained  from 
the  marginals  gets  closer  to  the  true  final  state. 

The  results  discussed  in  the  previous  section  were  obtained  as  steady 
states  of  this  modified  process,  using  a  time  window  of  10  iterations  (for 
the  experiment  of  Fig.  8,  an  additional  window  of  100  iterations  was 
used).  In  till  cases,  the  steady  state  was  obtained  after  approximately  450 
iterations.  Fig.  10  shows  snapshots,  taken  every  50  iterations,  of  the  state 
of  the  line  process  for  the  surface  reconstruction  problem  on  a  100  X 
100  lattice,  with  2000  observations,  portraying  a  central  square  figure  at 
constant  height  over  a  uniform  background. 


Fig.  8.  (a)  Observations  of  a  tilted  square  (slope  =  0.2)  over  a 
background  at  height  1.0  with  added  white  Gaussian  noise  (cr  = 
0.1).  White  pixels  denote  missing  observations,  (b)  Boundaries 
found  by  the  Algorithm,  (c)  Interpolation  with  no  boundaries, 
(d)  Reconstructed  surface. 
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Fig.  10.  Snapshots  of  the  Annealing  process. 


5.  Discussion. 


The  results  we  have  presented  indicate  the  plausibility  of  using  varia¬ 
tional  principles  that  include  our  prior  knowledge  about  the  geometry  of 
the  discontinuities  of  a  piecewise  smooth  surface,  to  reconstruct  it  from 
sparse  and  noisy  data,  preserving  the  boundaries  between  continuous 
subregions.  This  method  may  also  have  a  more  general  significance,  since 
most  problems  in  early  vision  are  ill-posed  and  must  be  regularized  [P2]. 
The  regularization  methods  used  so  far  are  based  on  smoothness  assump¬ 
tions,  and  for  this  reason,  all  have  problems  when  dealing  with  functions 
that  are  only  piecewise  continuous.  In  this  sense,  the  approach  we  have 
presented  may  be  regarded  as  an  extension  of  standard  regularization 
methods  for  handling  discontinuities  (see  [P2]).  The  algorithms  presented 
are  admittedly  slow.  Their  computational  efficiency,  however,  can  be 
greatly  improved  by  the  use  of  parallel  hardware,  particularly  in  view  of 
the  local  nature  of  the  support  for  the  state  updating  operations  (equations 
(9)  and  (10)).  It  would  be  interesting  to  investigate  the  implementation  of 
this  algorithm  on  the  "Connection  Machine",  currently  under  construction 
at  the  A.I.  laboratory  [HI]  (see  also  [FI]). 

5.1  Relation  to  other  Boundary  Preserving  Surface  Reconstruction  Techniques. 

The  use  of  MRF  models  is  closely  related  to  the  use  of  physical 
models  for  the  surfaces  that  one  wants  to  reconstruct.  Thus,  for  example, 
Terzoupulos  [T2]  proposes  the  use  of  a  thin  plate  model  to  embody  the 
knowledge  about  the  smoothness  of  the  surface.  A  threshold  on  the 
bending  moment  of  this  plate  is  then  used  to  locate  the  discontinuities, 
and  the  plate  is  allowed  to  break  at  these  points.  This  technique  is 
computationally  more  efficient,  since  the  energy  functionals  involved  are 
convex;  however  the  use  of  the  method  we  are  proposing  seem  to  have 
some  definite  advantages: 

(i)  From  a  conceptual  viewpoint,  it  is  better  to  perform  the  inter¬ 
polation  and  boundary  detection  tasks  at  the  same  time,  rather 
than  approximating  an  everywhere  smooth  surface  first,  since 
this  operation  hides  the  discontinuities  that  one  then  tries  to 
find  in  the  second  phase. 

(ii)  In  our  method,  the  values  of  the  parameters  depend  only  on 
the  average  height  of  the  jumps  that  one  wants  to  consider 
as  boundaries  in  the  reconstructed  surface ,  and  thus,  they  are 
independent  of  the  location  of  the  observations.  If  these  are 
sparsely  located,  the  bending  moment  of  the  thin  plate  may 
be  low,  even  when  the  discontinuity  is  relatively  large,  and  the 
threshold  method  may  fail. 


(iii)  A  priori  knowledge  about  the  shape,  orientation  and  position  of 
the  discontinuities  can  be  easily  incorporated  by  choice  of  the 
potentials  of  the  line  process. 

(iv)  The  same  algorithm  can  be  used  for  surface  interpolation,  noise 
elimination  (smoothing)  and  boundary  detection. 

5.2  Extensions  and  Open  Problems. 

This  method  can  be  easily  extended  to  the  case  of  piecewise  smooth 
(not  necessarily  straight)  boundaries  by  using  a  4-valued  line  process 
to  include  edge  orientation  (see  [Gl]).  This  extension  will  increase  its 
practical  value,  and  permit  its  use  in  a  variety  of  Computer  Vision  tasks 
(processing  of  depth,  motion  and  brightness  data)  as  well  as  in  other  fields 
(geophysics,  medicine,  etc.) 

Another  interesting  extension  is  the  use  of  the  knowledge  about  the 
position  of  the  boundaries  provided  by  some  process  (for  example,  the 
output  of  a  conventional  edge  detector  operating  on  the  components  of  a 
stereo  pair)  to  influence  the  reconstruction  of  a  different  related  surface 
(such  as  the  disparity  surface  obtained  from  a  stereo  matcher).  This  can 
be  easily  accomplished,  in  principle,  by  modifying  the  potentials  of  the 
line  process,  so  that  the  presence  of  a  related  edge  at  a  line  location 
lowers  the  energy  of  the  corresponding  configurations.  We  are  currently 
working  on  these  developments. 

An  important  problem  that  has  to  be  addressed  to  make  this  tech¬ 
nique  more  practical,  is  the  improvement  of  its  computational  efficiency. 
One  needs  to  accelerate  the  convergence  of  the  annealing  scheme,  and 
also,  see  if  it  is  possible  to  exploit  the  structure  of  the  energy  functional  of 
this  particular  problem,  to  develop  more  efficient  (possibly  deterministic) 
algorithms  to  find  its  global  minimum. 

Another  open  question  is  the  determination  of  the  optimal  values  for 
the  parameters  of  the  energy  functionals.  It  will  be  interesting  to  explore 
in  this  connection  the  use  of  statistical  methods  [B1,K3],  regularization 
techniques  [P2,T1,W1],  and  learning  algorithms  [H2J. 
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